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INTRODUCTION 



A. BACKGROUND 

In analyzing ocean acoustic pulse-propagation problems, 
accurate phase calculations must be performed in order to 
predict the acoustic signal at the receiver. Linear systems 
theory provides analytical expressions for analyzing the 
propagating acoustic field [Refs. 1-4]. These well-known 
equations form the basis for the FORTRAN programs developed in 
this thesis. 

Program input is depth versus sound-speed data pairs where 
the speed of sound is an arbitrary function of depth. The 
ocean is viewed as a linear, time-variant, space-variant 
filter. The WKB approximation can specify this filter's 
transfer function [Refs. 1 - 4]. For the arbitrary sound- 
speed profile, no exact transfer function exists. The 
transfer function requires a method for evaluating signal 
phase at the receiver. 

Two phase evaluation methods are presented and contrasted. 
Each method calculates the phase of the acoustic signal for a 
specified horizontal range. The first method overlays a 
piecewise linear profile on the input sample values for the 
speed of sound at various depths. The theory of ray acoustics 
is used to calculate travel time and phase. The second method 
fits a smooth cubic spline curve to the input samples. This 
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method solves the propagation problem using a system of three, 
first-order differential equations [Ref. 5]. Both methods 
make phase calculations using signal travel times. 
Additionally, each method is capable of identifying eigenrays 
or rays that directly connect the signal source to the 
receiver. 

Another method is introduced to validate results. Direct 
numerical integration is performed to calculate phase for a 
single gradient, linear, sound-speed profile. A separate 
FORTRAN program implements this numerical integration routine 
providing a totally independent verification. 

An overview of the theory behind each method is presented 
in Chapter II. Computer simulation results are presented in 
Chapter III for various input sound-speed profiles. Tabular 
results of phase calculations are presented for each of the 
three analysis methods. Using input from the piecewise linear 
and the cubic spline/differential equation solution methods, 
ray traces of the propagating field are presented as a visual 
aid in interpreting the results. 
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II. THEORETICAL BACKGROUND FOR THE 
EVALUATION OF PHASE INTEGRALS 

A. UNDERWATER ACOUSTIC PULSE PROPAGATION 

Ocean acoustic pulse-propagation models can be derived by 
using the principles of linear, time-variant, space-variant, 
systems theory and the physics of wave propagation in 
inhomogeneous media. Linear systems theory allows for the 
development of an ocean medium transfer function. An ocean 
medium transfer function that is based on the WKB 
approximation has been derived and is given by References 1 
and 2 as follows: 



H(f, f r ,y 0 ;y) 



A 

/rwr 



-J f MOdC 
J»o 



e J 2 »rf,(y-y 0 ) 



( 2 . 1 ) 



where 



f is the frequency in Hz 

f r is the radial, spatial frequency in 
cycles per meter 

y 0 is the source depth in meters 



A 



-d 



K (y 0 ) 



4 7T f,, 



( 2 . 2 ) 
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ky(y) is the propagation vector component in 
the Y direction with units of radians 
per meter and is given by 



My) 



'± 2n ( [f/c(y) ] 2 - f r 2 }* ,f r < f/c (y) (2.3) 

* j27r(f r 2 - [f/c(y)] 2 } % ,f r > f/c (y) (2.4) 



c(y) is the sound-speed expressed as a function 
of depth with units of meters per second 

f is the transmitted (input) spatial frequency 

in the Y direction at the source in cycles per 
meter and is given by 



c 0 is the speed of sound in meters per second at 
the source depth y 0 , that is, c 0 = c(y 0 ) ; 

The plus (minus) sign in Eqs. (2.3) and (2.5) is chosen 
whenever y - y 0 > 0 (y - y 0 < 0) . The minus (plus) sign in 
Eqs. (2.4) and (2.6) corresponds to the plus (minus) sign in 
Eqs. (2.3) and (2.5). This thesis deals only with propagating 
waves. Therefore, Eqs. (2.4) and (2.6) representing the 
generation of evanescent waves (i.e., decaying exponentials), 
will not be used. 

The evaluation of the phase integral in Eq. (2.1), namely, 



± [ (f/c 0 ) 2 - f r V ! , f r < f/c 0 

T [f r 2 - (f/c 0 ) 2 ]* , f r > f/c 0 



(2.5) 



( 2 . 6 ) 



and 




(2.7) 
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has been attempted in past studies in which solution 
techniques included direct integration and binomial expansions 
[Ref. 3]. 

Computer programs using these techniques proved to be very 
expensive in terms of CPU time. This thesis compares four 
different methods for evaluating the phase integral. The four 
methods are presented in the following sections. 

B. DIRECT INTEGRATION 

Computer code was developed to directly evaluate the phase 
integral given by (Eq. 2.7) using numerical integration 
routines from the International Math and Statistics Library 
(IMSL) . Direct integration is a time-consuming technique that 
was applied only under the following constraints: 

- free-space acoustic propagation 

- linear sound-speed profile with a single gradient g 
in inverse seconds and 

- propagating signals that have not passed through a 
turning point in the medium. 

Expressions for the propagation vector component ky(y) of 
the phase integral are given by Eqs. (2.3) and (2.4). For a 
linear sound-speed profile with a single gradient, the 
dependence of c(y) on depth y is given by 

c (y) = c(y 0 ) + g (y - y 0 ) (2.8) 

where 

c(y 0 ) is the sound-speed at the source depth 
y 0 meters 
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g is the sound-speed gradient, and 

y is the desired depth for sound-speed. 

Since the direct integration approach is limited to free- 
space, single gradient, propagation problems, it is used only 
to validate the results of other solution techniques for 
several simple test cases. The results will be compared 
against the next two methods to be presented. These methods 
will apply the theory of ray acoustics [Ref. 4 & 5] to find 
travel time. Phase calculations will be based on travel time 
calculations . 

C. PHASE CALCULATIONS BASED ON PIECEWISE LINEAR SOUND-SPEED 

PROFILES 

A FORTRAN program developed by Lim [Ref. 5] was the ray 
acoustics propagation code that was used to draw ray trace 
plots and to calculate travel time, in addition to the depth, 
angle of propagation and path length along a ray path as a 
function of horizontal range. This method applies ray 
acoustics to a piecewise linear sound-speed model of the ocean 
medium. Based on the input sound-speed versus depth data 
pairs, layers are defined in the ocean medium. Each layer has 
an upper and lower boundary at specific depths. A constant 
sound-speed gradient g is calculated for each layer. The 
sound-speed for any desired depth is computed using Eq. (2.8) 
with the appropriate gradient. 

In a constant gradient medium, ray acoustics theory allows 
calculation of travel time with closed-form equations. For an 
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incremental increase in horizontal range from the source, the 
angle of arrival for the propagating ray can be shown to be 
[Ref. 5], 

j8(y) = cos 1 [cos/3 0 - b*g(y) * rngstp] (2.9) 

where 

/3 0 is the ray launch angle, 

g(y) is the gradient in the layer at depth y, and 

rngstp is the incremental increase in horizontal 
range. 

The ray parameter b in Eq. (2.9) is given by [Refs. 4 & 5] . 

b = sin 0 O / c(y 0 ) (2.10) 

Knowing the arrival angle /3(y), the ray depth y is given by 
[Refs. 4 & 5] 



y 




sinff(y) 
sin/3 (y 0 ) 



( 2 . 11 ) 



The arrival angle also allows calculation of travel time from 
[Refs. 4 & 5] 



T 



l ln tan[/3(y)/2] 
g tan|/3 ( y 0 ) / 2 J 



( 2 . 12 ) 



Once the travel time (in seconds) is found for the desired 
horizontal range, total phase is given by 

0 T = 27TfT (2.13) 

The phase integral given by Eq. (2.7) represents a phase 

change in the depth (y) direction. The constant value of the 
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propagation vector component in the radial direction allows 
for easy calculation of phase change in the radial direction. 
This value is given by 

0 R = 27rf r * hrzrng (2.14) 

where hrzrng is the total horizontal range traveled. 

The phase in the depth direction is 

e y = e T - e R . (2.15) 

Phase expressed in radians is a modulo 2n function. The 
proper solution for the phase integral Eq. (2.15) must be 
expressed as a modulo 2ir function: 

0 y = modulo (0 y , 27r ) . (2.16) 

This method (referred to as method 1 in Ref. 5) is capable 
of analyzing propagating rays at any horizontal range in 
either free space or a bounded medium. 

D. PHASE CALCULATIONS BASED ON AKIMA CUBIC SPLINES AND 
ORDINARY DIFFERENTIAL EQUATIONS 

This application of ray acoustics applies a more 
sophisticated treatment to the depth versus sound-speed data 
pairs sampled from the ocean medium. The data pairs are used 
to form a smooth Akima cubic spline sound-speed profile. 
Splines offer the advantage of using all data points in 
generating a profile, and they place no restrictions on 
spacing between data points. The Akima version of the cubic 
spline was chosen for its excellent ability to combat wiggles 
in the profile, that is, it suppresses oscillations that would 
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cause overshoots and undershoots in the sound-speed versus 
depth profile. 

This method uses ray acoustics theory to generate a system 
of three, first-order differential equations [Ref. 5]: 

Yi = y 2 (2.17) 



-c (y^ 

b 2 c 3 (y : ) 



(2.18) 



y 3 = = (2.19) 

i bc 2 ( 7l ) 

where 

c(y 1 ) is the derivative of the sound-speed with 
respect to depth at depth y 1 , 

y 1 is the ray depth, 

y 2 is the cotangent of the ray's arrival angle 
0(y), and 

y 3 is the travel time of the ray. 

Once solved, this system of differential equations allows 
phase to be calculated from travel time. The phase 
calculation is performed exactly as shown for the piecewise 
linear sound-speed profile. 

Like the previous method, the Akima cubic 
spline/differential equation solution can be used on free 
space and bounded media problems. 
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E. EIGEN RAYS 

Eigenrays are propagating rays that exactly connect the 
sound source to the receiver. The FORTRAN propagation code 
developed can search for and identify eigenrays. The input 
required is 

- depth versus sound-speed data pairs, 

- source depth in meters, 

- receiver depth and range in meters, 

- angle step in degrees between possible eigenrays to 
be evaluated, and 

- allowed depth error denoted y error in meters. 

Rays passing within the allowed error or tolerance y error of the 
receiver are identified as eigenrays. The eigenray mode can 
employ either the piecewise linear sound-speed profile or the 
Akima cubic spline/differential equation method for ray 
propagation. 

The eigenrays are found by trial-and-error . This search 
method was chosen after the IMSL DBVPMS program failed to 
solve the problem. The DBVPMS program is a differential 
equation solver that was applied to the system of differential 
equations given by Eqs. (2.17) through (2.19). The routine 
uses the shooting method to find eigenray solutions to 
boundary value problems. It was unable to converge to a 
solution. Solutions for the acoustic problems investigated 
are difficult for this algorithm because of 
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- long propagation distances, 

- the inhomogeneous ocean presents a continuously 
varying medium with discrete boundaries, and 

- the system of differential equations has non-constant 
coefficients whose values are a function of depth. 

When the eigenray mode is selected, rays that are not 

refracted or reflected to within y error of the receiver are 

discarded . 
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Ill 



COMPUTER SIMULATION RESULTS 



A. OVERVIEW 

The computer simulation results presented in this 
chapter perform phase calculations intended to: 

- show that ray acoustics can be used to evaluate 
the phase integral producing the same results as 
direct numerical integration, 

- show that the two methods of applying ray acoustics 
theory (presented in Chapter II) produce equivalent 
results within an acceptable tolerance, 

- compare speed versus accuracy for the two methods for 
a variety of test cases, and 

- demonstrate the ability to identify eigenrays using 
either of the ray acoustics methods. 



B. DIRECT INTEGRATION 

In the previous chapter, 0 y (y) was shown as the solution 
to the phase integral in the depth direction. 



e ' (y) = iIv C)dC 



( 2 . 7 ) 



The straight-forward evaluation technique is the brute 
force approach using direct numerical integration. The phase 
integral will be evaluated using direct numerical integration 
for a simple test case. The results will then be compared to 
phase values obtained using ray acoustics calculating phase 
values from travel times. 
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Conditions for the test case are: 



- source depth of 10 meters, 

- receiver depth of 100 meters, 

- linear sound-speed prpfile with a single, constant 
gradient of 0.016 sec and 

- four values of radial spatial frequency (FR values) 
will be evaluated. 

For each FR value, tabular results will show the 
corresponding ray launch angle /? ( y 0 ) . The relation between 
the FR value and this launch angle is [Ref. 4]: 

0(y o ) = sin 1 C FR * c(y 0 )/f] (3.1) 

where 

f is the frequency in HZ, and 

c(y 0 ) is the speed of sound at the source depth. 
Table 3 . 1 shows the results produced by direct numerical 
integration. Tables 3.2 and 3.3 show the results for the 
Akima cubic spline/ordinary differential equation solver 
(ODE solver) and the piecewise linear approximation ray 
acoustics techniques, respectively. 

The parameter of interest is THETAY representing the phase 
change in the depth or Y direction. Table 3.1 shows these 0 
values to be approximately 92, 86, 75, and 56 radians for the 
four FR values chosen. Tables 3.2 and 3.3 contain four sets 
of entries, each corresponding to an FR value. These entries 
begin at 0 meters range and end at the range corresponding to 
100 meters depth; thus, solving the phase integral for the 0 y 
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TABLE 3.1 

PHASE INTEGRAL SOLUTION BY DIRECT NUMERICAL INTEGRATION 



INPUT DATA FOR PHAGE INTEGRAL EVALUATION 



F - 250.0 HZ 

EA0S - 0.000000 
YREF • 0.0 M 



YL ■ 10.0 M YU ■ 100.0 M RATIO ■ 0.9999 NTR « 5 

EREL ■ 0.000100 

CYREF - 1500.0 M/SEC G - 0.016000 (1/SEC) 



METHOD ■ 2 



EVALUATION BY 1MSL10 ROUTINE DODAG WITH 1RULE • 2 



FR ( CYCLES/M ) BETAYO(DEG) RANGE AT 

* YU 



ANGLE OF ARRIVAL 
BETAY ( DEG ) 



TURNING POINT TURNING POINT TMETAR(RAD) THETAY(RAD) 
DEPTM(M) RANGE (M ) 



0.0332944859 


11.525 


18.36009 


11.536 


0.0665889718 


23.552 


39.25213 


. 23.576 


0.0998834576 


36.824 


67.43892 


36.866 


0.1331779435 


53.049 


119.80722 


53.122 



0 . 37555E*06 
0.14090E*06 
0. 62682E*05 
0.23S74E*05 



0.45984E*06 
0 . 2 15 1 OE* 0 6 
0. 12522E*06 
0. 7052?E*05 



0.38408E*01 
0 . 1 642 3E *02 
0 . 42 3 24E*02 
0. 10025E*03 



0.92291679S6E*02 
0.8633840249E*02 
0.75378S7l91E*02 
0 . 565 7 3684 1 3E* 02 



TABLE 3.2 

PHASE INTEGRAL SOLUTION BY ODE SOLVER 



ERAYS « F 
YO • 10.0 M 

FREOC/CMAX « 0.16649 CYCLES/M 
DLTFR » 0.033294 CYCLES/M 

NOTE 1 FRM1N » DLTFR 



PRTALL = T 
DEPTH « 100.0 M 

RATIO « 0.9999 
RNGSTP • 1 .0 M 



FREOC • 250.0 HZ CMAX « 1501.6 M/SEC 

FRMAX - 0.16647 CYCLES/M NFR$ « 5 

MRZRNG ■ 150.0 M YR » 40.0 M 



FR - 0. 332944058B1726E-O1 CYCLES/M 

BET AO « 1 1.5245752 1170979 DEG 



RANGE (M) 


DEPTH(H) 


TRVLT (SEC ) 


BETAY(DEG) 


BTEST 


THETAT(RAD) 


THETAR (RAD) 


THETAY(RAD) 


0.0 


10.00 


0.000000 


1 1.525 


0. 13318E-03 


0. 0000000E*00 


O.OOOOOOOE'OO 


0.0000000E*00 


18.4 


100.00 


0.061200 


11.536 


0.133 18E-03 


0 . 96 1 3253E* 02 


0.3840846E*01 


0 . 92291 68E*02 



FR * 0.665889717634S2E-01 CYCLES/M 

BETAO » 23.55170283980132 DEG 



RANGE (M ) 


DEPTH(M) 


TRVLT ( SEC ) 


BETAY(DEG) 


BTEST 


THETAT(RAD) 


THETAR (RAD) 


TMETAY(RAD) 


0.0 


10.00 


0.000000 


23.552 


0.26636E-03 


0.0000000E*00 


o.oooooooe*oo 


0 .OOOOOOOE‘00 


39 . 3 


100.00 


0.065420 


23.576 


0 .26636E-03 


0.1 02 761 IE *03 


0. 1642273E*02 


0 .8633840E*02 



FR * 0.998834576451 78E-01 CYCLES/M 

BETAO * 36.82440911644773 DEG 



RANGE (M) 


DEPTH(M) 


TRVLT ( SEC ) 


BET AY (DEG) 


BTEST 


THETAT (RAD ) 


THETAR ( RAD ) 


THETAY(RAD) 


0.0 


10.00 


0.000000 


36 .824 


0.39953E-03 


o ,oooooooe»oo 


0 .OOOOOOOE'OO 


o.oooooooe*oo 


67.4 


100.00 


0.074932 


36.866 


0. 39953 E- 03 


0.1 177023E*03 


0 . 42323 74E* 02 


0 , 753 785 7E* 02 



FR * 0 . 1331 7794 552690E*00 CYCLES/M 



BETAO = 


53.0492 


8575821315 DEG 












RANGE(M) 


DEPTH(H) 


TRVLT (SEC) 


BETAY(DEG) 


BTEST 


THETAT (RAD ) 


THETAR (RAD) 


THETAY (RAD ) 


0.0 


10.00 


0.000000 


53.049 


0 .5327 1 E-03 


o . oooooooe*oo 


O.OOOOOOOE'OO 


o . oooooooe*oo 


119.8 


100.00 


0.099839 


53.122 


0 .53271E-03 


0. 1568262E*03 


0. 1 002525E * 0 3 


0.56573 68E*02 
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TABLE 3.3 

PHASE INTEGRAL SOLUTION BY PIECEWISE LINEAR APPROXIMATION 



RAY TRACING USING PIECEWISE LINEAR APPROXIMATION 
CASE 1 4 LA 



NDATA = A 

G( 1 ) = 0 . 16000 E- 01 1/SEC 

GC 3 ) = 0 . 1 6000E-0 1 1/SEC 



NUMBER OF GRADIENTS = 3 

GC 2 ) * 0. 16000 E- 01 1/SEC 



ERAYS = F 
Y0 = 10 .0 M 

FREOC/CMAX = 0.166A9 CYCLE: 
DLTFR = 0 . 0 3 32 9 A CYCLES/M 
NOTE: FRMIN = DLTFR 



PRTALL = T 
DEPTH = 100.0 M 

/M RATIO * 0.9999 
RNGSTP * 1.0 M 



FREOC «= 250.0 HZ 

FRMAX = 0 . 1 6 6 A 7 CYCLES, 
HR2RNG * 150.0 M 



CMAX ■ 1501.6 M/SEC 

NFRS * 5 

YR * A 0 . 0 M 



FR * 0. 3329AA85881726E-0 1 CYCLES/M 

BETA0 ■ 1 1 .52A57521 170980 DEG 



RANGE CM) 


DEPTH ( M ) 


TRVLT (SEC ) 


BETAY ( DEG ) 


B TEST 


THETAT(RAD) 


THETAR(RAD) 


THETAY(RAD) 


0.0 


10.00 


0.000000 


11.525 


0 . 13318E-03 


0.0000000E 4 C0 


0 . 0000000E 4 00 


0 . 0000000E 4 00 


18. A 


100.00 


0.061200 


11 .536 


0.13318E-03 


0 . 9 6 1 325 3E 4 02 


0 . 38 A 08 A 6E 4 0 1 


0. 92291 68E*02 



FR = 0. 


66588971763A52 


E-01 CYCLES/M 












EETAO = 


23.5517028 


3 98 0 1 3 A DEG 












RANGE (Ml 


DEPTH CM) 


TRVLT (SEC ) 


BETAY (DEG ) 


BTEST 


THETAT(RAD) 


THETAR(RAD) 


THETAYCRAD) 


0 . 0 


10.00 


0.000000 


23.552 


0.26636E-03 


0.0000000E 4 00 


0.0000000E 4 00 


0.OOOOOOOE 4 OO 


39.3 


100.00 


0.065A20 


23.576 


0.26636E-03 


0.102761 1E 4 03 


0 . 16A2273E 4 02 


0.863 38 A0E 4 02 



FR = 0.99883A576A5178E-01 CYCLES/M 



BETAO = 


36 .82AA091 1 6 A A 7 7 6 DEG 












RANGE(M) 


DEPTH (M ) 


TRVLT (SEC ) 


BET AY ( DEG ) 


BTEST 


THETAT(RAD) 


THETAR (RAD) 


THETAY(RAD) 


0.0 


10.00 


0 .000000 


36 . 8 2A 


0. 39953 E- 03 


0 .0000000E*00 


0 . 0000000E 4 00 


0.0000000E 4 00 


67. A 


100.00 


0 . 0 7 A 9 32 


36.866 


0 . 39953E-03 


0 . 1 1 77023E 4 03 


0.A23237AE*02 


0.753 785 7E*02 



FR » 0.1 


331779A352690E 4 00 CYCLES/M 












BETAO = 


53 . 0A92857582 1 320 DEG 












RANGE (M ) 


DEPTH (M ) TRVLT (SEC ) 


BETAY ( DEG ) 


BTEST 


THETATCRAD) 


THETAR (RAD) 


THETAY(RAD) 


0.0 


10.00 0.000000 


53.0A9 


0 .53271E-03 


0.0000000E 4 00 


0 .0000000E 4 00 


0.0000000E 4 00 


119.8 


100.00 0.099839 


53 . 122 


0. 53271 E- 03 


0 . 15682 62E 4 03 


0 . 1 0 02525E 4 0 3 


0. 56573 68E 4 02 
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value when a ray travels from 10 to 100 meters depth. The 
significance of these three tables is that the ray acoustics 
methods are yielding the same phase values 0 y as the direct 
integration method. 

The conclusion to be derived from this is that ray 
acoustics theory can be used to accurately evaluate the phase 
integral. 

C. VARYING SOUND-SPEED GRADIENTS 

Both ray acoustics methods will now be applied to three 

different test cases. A zero gradient, positive gradient 

-1 . -1 
(+0.016 sec ), and a negative gradient (-0.016 sec ) linear 

sound-speed profile will be used. These simulations will show 

that both methods produce equivalent results within an 

acceptable tolerance for a variety of media. 

The conditions common to each simulation run are: 

- source depth of 10 meters, 

- horizontal range of 1 kilometer, 

- ocean depth of 100 meters, 

- speed of sound at the surface of 1500 meters per 
second, and 

- five values of radial spatial frequency (FR) 
evaluated. 

The results are presented in Tables 3.4 through 3.9. The 
phase integral solution is the phase change in the Y (depth) 
direction listed as the MODULO TWOPI THETAY value in radians. 
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TABLE 3.4 

ZERO GRADIENT , PIECEWISE LINEAR SOLUTION 



RAY TRACING USING PIECEWISE LINEAR APPROXIMATION 
CASE 1A 

NDATA - 4 NUM8ER OF GRADIENTS - 3 

G( 11 • 0 . 0 00 00E* 00 1/SEC G( 2) - 0.00000E*00 1/SEC 

G ( 31 - 0 . 00000E+00 1/SEC 

Y0 • 10.0 M DEPTH - 100. 0 M FREOC ■ 250.0 HZ CMAX ■ 1500.0 H/SEC 

FREOC/CMAX ■ 0.U667 CYCLES/M RATIO ■ 0.9999 FRMAX « 0.16665 CYCLES/M NFRS • 5 

DLTFR > 0.033330 CYCLES/M RNGSTP • 1.0 M HRZRNO ■ 1000. 0 M 

NOTE: FRM1N > DLTFR 



NFR 


BETAO(DEG) 


DEPTM(M) 


TRVLT (SEC) 


8ETAY(DEG1 


THETAT(RAD 1 


THETAR (RADI 


THETAY(R AD 1 


MODULO TWOP1 
THETAY(RAD) 


l 


11.536 


90.51 


3.333667 


168.464 


0 .523651 1E 4 04 


0 . 2094 186E*03 


0 .5027093E 4 04 


0.5445951 


2 


23.576 


98.64 


1.666833 


156.424 


0.26 18256E*04 


0 . 4 188 3 7 IE*0 3 


0.2199419E*04 


0.3037135 


3 


36.866 


56.46 


1.111222 


143.134 


0.1745504E*04 


0.6282557E*03 


0.111 7248E*04 


5.1243041 


4 


53.122 


39.79 


0.833417 


126.878 


0 . 1 309 128E*04 


0.8376743E*03 


0.47 14536E 4 03 


0.2146886 


5 


89.190 


24.14 


0.666733 


89.190 


0. 1047302E»04 


0. 1047093E*04 


0 . 2094500E*00 


0.2094500 


TOTAL CPU TIME 


• 0 MIN 


, 22.63 SEC 













TABLE 3.5 

ZERO GRADIENT , ODE SOLVER SOLUTION 



RAY TRACING USING AKlMA CU81C SPLINE A ODE SOLVER 
CASE 2 A 
NDATA • 4 



Y0 ■ 10.0 M 

FREOC/CMAX • 0.16667 CYCLES/M 
DLTFR - 0.033330 CYCLES/M 

NOTE: FRM1N - DLTFR 



DEPTH • 100.0 M 

RATIO • 0.9*99 
RNGSTP • 1 .0 M 



FREOC • 250.0 HZ CMAX 
FRMAX - 0.16665 CYCLES/M NFRS 
HRZRNG • 1000.0 M 



1500.0 M/SEC 
5 



NFR 


8ETA0 (DEG ) 


DEPTH(M) 


TRVLT (SEC) 


8ETAY(DEG ) 


THETAT(RAD) 


THETAR (RAD) 


THETAY(RAD) 


MODULO TWOPI 
THETAY(RAD) 


1 


11.536 


90.51 


3.333667 


168.464 


0.5236511E*04 


0 .2094186E*03 


0 . 502709 3E*04 


0 . 5445951 


2 


23.576 


98.44 


1 . 666833 


156.424 


0.26 18256E*04 


0 . 4 188 3 7 1 E* 0 3 


0 . 21 994 19E*04 


0.3037135 


3 


36.866 


56.46 


1.111222 


143.134 


0. 1745504E*04 


0 . 6282557E*03 


0.1 11 7248E*04 


5.1243041 


4 


53.122 


39.79 


0.833417 


126.878 


0. 1309128E*04 


0.8376743E*03 


0.471 4536E*03 


0.2146886 


5 


89.190 


24.14 


0.666733 


89.190 


0. 1047302E*04 


0. 1047093E*04 


0.2094500E*00 


0. 2094500 



TOTAL CPU TIME • 15 MIN , 52.38 SEC 
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TABLE 3.6 

POSITIVE GRADIENT, PIECEWISE LINEAR SOLUTION 



RAY TRACING USING PIECEWISE LINEAR APPROX lHAT ION 
CASE 18 

NDATA - 4 NUH8ER OF GRADIENTS ■ 3 



G ( 1 ) 


* 0 . 1 4000E-0 1 


l/SEC G ( 


2 > ■ 0. 1 4 OOOE-O 1 


l/SEC 






G ( 3) 


* 0 . 1 4000E-0 1 


1/SEC 










Y0 > 


10.0 M 




DEPTH ■ 100.0 M 


FREOC ■ 250.0 H2 


CHAX » 


1501.4 M/SEC 


FREOC/CMAX • 0.1664« 


CYCLES/H 


RATIO * 0.9999 


FRHAX - 0.14447 CYCLES/M 


NFRS • 


5 


DLTFR 


- 0.033294 CYCLES/M 


RNGSTP • 1 .0 H 


MRZRNG • 1000.0 H 







NOTE: FRH1N ■ DLTFR 



NFR 


BETAO(DEG) 


DEPTM(M) 


TRVLT(SEC) 


8E T AY( DEG ) 


THETAT (RAD ) 


THETAR (RAD ) 


TMETAY(RAD) 


MODULO TWOPI 
THETAYtRADJ 


1 


11.525 


87.80 


3.333452 


148 .444 


0 . 52 34408E *04 


0.2091 954E *03 


0 .5027293E*04 


0.7442549 


2 


23.552 


94.99 


1 .444824 


154.425 


0.2418241E*04 


0.4 18 3 908E * 0 3 


0 . 2 1 9 9850E * 04 


0 .7348999 


3 


. 34.024 


55 • 38 


1.111192 


143.155 


0. 1745454E*04 


0. 427584 3E *03 


0. 1 1 1 7870E*04 


5.7442108 


4 


53.049 


38 .73 


0.833378 


124.927 


0. 1 309044E *04 


0.834781 7E*03 


0 . 4 7228 4 7E • 0 3 


1.0457792 


5 


87.343 


50.7 1 


0.447000 


87.974 


0 . 1047721E»04 


0. 1045977E*04 


0. 17437 98E*01 


1 .7437978 


TOTAL CPU TIME 


• 0 HIM 


. 22.05 SEC 













TABLE 3.7 

POSITIVE GRADIENT, ODE SOLVER SOLUTION 



RAY TRACING USING AKlHA CUBIC SPLINE l ODE SOLVER 
CASE 28 
NDATA ■ 4 

YO • 10.0 H DEPTH • 100.0 H FREOC ■ 250.0 M2 CHAX ■ 1501.4 M/SEC 

FREOC/CMAX * 0 . 1 6 6 <S « CYCLES/M RATIO » 0.9999 FRHAX ■ 0.1444? CYCLES/M NFRS - 5 

DLTFR • 0.033294 CYCLES/H RNGSTP • 1.0 H HR2RNG • 1000.0 H 

NOTE: FRM1N • DLTFR 



NFR 


BETAO(DEG) 


DEPTH! H> 


TRVLT(SEC) 


SET AY! DEG > 


THETAT ( RAD ) 


THETAR(RAD) 


THETAYCRAD) 


MODULO TWOP 
THETAY(RAD) 


1 


1 1.525 


87.80 


3.333452 


148 .444 


0 .52 34408E *04 


0.2091 954E * 03 


0.502?293E*04 


0.7442493 


2 


23.552 


94.99 


1.444824 


154.425 


0.2418241E*04 


0 .4183908E*03 


0 .2 1 99850E *04 


0.7349005 


3 


34.824 


55.38 


1.111192 


143.155 


0. 1745454E*04 


0. 4275843E*03 


0. 11 I 78 70E*04 


5.7442125 


4 


53.049 


58.73 


0.83:378 


124 . 927 


0. 1309044E*04 


0.834781?E*03 


0.4 7228 47E*03 


1.0457793 


5 


87.343 


50.71 


0.447000 


87.974 


0. 1047721E»04 


0.1045977E*04 


0.1 743798E*01 


1.7437978 



TOTAL CPU TIME • 14 MIN . 24.30 SEC 
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TABLE 3.8 

NEGATIVE GRADIENT, PIECEWISE LINEAR SOLUTION 



RAY TRACING USING PIECEWISE LINEAR APPROXIMATION 
CASE 1C 

NDATA • A NUM8ER OF GRADIENTS • 5 

G( 1) « -0.1 6000E- 0 1 1/SEC G( 2) « -0.16000E-01 1/SEC 
G( 5) » -0.1 6000E-0 1 1/SEC 



YO * 1 0 . 0 M 

FREOC/CMAX • 0.16667 CYCLES/M 
DLTFR • 0.053530 CYCLES/M 

NOTE: FRM1N ■ DLTFR 



DEPTH - 100.0 M 

RATIO * 0.9999 
RNG3TP » 1 .0 M 



FREOC - 250.0 HZ CMAX . 1500.0 M/SEC 

FRMAX « 0.16665 CYCLES/M NFRS - 5 
HRZRNG • 1000.0 M 



NFR 8ETA0 (DEG) DEPTM(M) TRVLT(SEC) 8ETAYCDEG) THETAT(RAD) 



1 


11.535 


87.78 


3.337240 


168.475 


0 . 5242 1 24 E* 04 


2 


23.573 


96.98 


1 .668622 


156.450 


0 . 262 1 065E *04 


3 


36.861 


55.32 


1 . 1 1 2 A 58 


143.160 


0.1 747414E*04 


A 


5 3 . 1 1 A 


38 .63 


0.834J45 


126.909 


0. 1310586E*04 


5 


88.835 


35.67 


0.667046 


88.224 


0. 1047793E*04 


TOTAL CPU TIME 


» 0 MIN 


. 22.61 SEC 







THETAR (RAD ) 

0.209A186E»03 
0.4 188371E*03 
0 . 6282557E*03 
0 .8 5 767 43E*03 
0 . 1047093E*04 



THETAY(RAD) 

0.5032706E*04 
0.2202228E*04 
0 . 1 1 1 9 1 59E* 04 
0 . 4 729 1 1 7E* 03 
0 . 6 9996 12E* 00 



MODULD TWOP1 
THETAY(RAD) 
6.1576516 
5 . 1 1 34 65 3 
0.7515713 
1.6727921 
0.6999612 



TABLE 3.9 

NEGATIVE GRADIENT, ODE SOLVER SOLUTION 



RAY TRACING USING AK1MA CU81C SPLINE & ODE SOLVER 
CASE 2C 
NDATA * 4 

YO * 10.0 M DEPTH * 100. 0 M FREOC « 250.0 HZ CMAX * 1500.0 M/SEC 

FREOC/CMAX • 0.16667 CYCLES/M RAT1D * 0.9999 FRMAX » 0.16665 CYCLES/M NFRS * 5 



DLTFR 


* 0. 033330 CYCLES/M RNGSTP * 1. 


0 M HRZRNG 


• 1000.0 M 






NOTE: 


FRM1N « 


DLTFR 














NFR 8ETA0(DEG) 


DEPTH(M) 


TRVLT(SEC) 


8ETAYCDEG) 


THETAT(RAD) 


THETAR (RAD) 


THETAY(RAD) 


MDDULD TWOP1 


















THETAY(RAD) 


1 


11.535 


87.77 


3.337240 


168.475 


0 . 5242 1 24 E ♦ 04 


0 .2094 186E*03 


0.5032706E*04 


6.1576656 


2 


23.573 


96.98 


1 .6 68 622 


156.450 


0.2621 065E*04 


0.4 188 371E*03 


0 . 2202228E* 04 


3.1134683 


3 


36 .86 1 


55.32 


1 . 1 124 38 


143.160 


0. 1747414E»04 


0.6282557E*03 


0.1119 159E*04 


0.7515729 


4 


53.114 


38.63 


0 . 8 3 4 3 4 S 


126.909 


0. 1310586E*04 


0. 837674 3E *03 


0 .4729 1 1 7E*03 


1.6727921 


5 


88 .835 


35.67 


0.667046 


88.224 


0. 1047793E*04 


0 . 1047093E*04 


0. 699961 2E*O0 


0.6999612 



TOTAL CPU TIME - 16 MIN , 21.26 SEC 
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Tables 3.4 and 3.5 show the results for the zero gradient 
case. A comparison between two tables shows that the phase 
values (MODULO TWOPI THETAY) agree perfectly. Additionally, 
all other calculated values in the tables agree. The only 
difference is the TOTAL CPU TIME. Both methods are using a 
conservative range step (RNGSTP) of 1 meter, i.e., the 
incremental increase in range is 1 meter in the propagation 
calculations. For these conditions, the ODE solver is taking 
over 30 times longer to run the simulation. 

Tables 3.6 and 3.7 show the results for the positive 
gradient case. A comparison between these two tables, as 
well as the negative gradient results shown in Tables 3.8 and 
3.9, shows acceptable agreement between the two very 
different calculation methods. The only discrepancies in 
calculated values occurs in the modulo 2 n 0 y phase values 
needed for solving the phase integral. The phase calculation 
is most challenging because of the nature of phase (being a 
modulo 27 r function) . The discrepancies for both the positive 
and negative gradients are on the order of hundred-thousandths 
of a radian or less. The CPU times continue to follow the 
pattern seen earlier that the ODE solver requires over 3 0 
times longer to complete a simulation run. 

The following ray trace plots, Figures 3.1 through 3.6, 
correspond to Tables 3.4 through 3.9. The plots assist 
interpretation of the data, but as expected from the tabular 
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RAY TRRflF. : PIECEWISE LINERR RPPROXIMRTIQN 




Figure 3.1 Ray trace corresponding to Table 3.4 
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HRZRNG : 1.0 KM 




Figure 3.2 Ray trace correseponding to Table 3.5 
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RflY TRACE: PIECEWISE LINEAR RPPROXIffflTmM 
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Figure 3.3 



Ray trace corresponding to Table 3.6 
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HRZRNG : 1.0 KM 



RRT TRACE: RKIMR CUBIC SPLINE j QDC SOLVER 
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Figure 3.4 Ray trace corresponding to Table 3.7 
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HRZRNG: 1.0 KM 



RAY TRACE : PIECEWISE LINEAR APPROXIMATION 




(W)Hld30 



Figure 3 . 5 



Ray trace corresponding to Table 3.8 
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HRZRNG : 1.0 KM 



RflY TRRCEj RKW CUBIC 




Figure 3 . 6 



Ray trace 



corresponding to Table 



3.9 



26 



HRZRNG: 1.0 KM 



results, no variations in depth or angle of arrival can be 
detected. 

D. CHALLENGING TEST CASES 

A smooth sound-speed profile in the form of a half 
sinusoid describes the medium for the next set of 
simulations. To this point, the agreement in accuracy between 
the two methods has been excellent; however, the shape of a 
sine curve is difficult to approximate with piecewise linear 
segments. While this will be a challenging test, smoothly 
curved profiles do occur naturally. The following set of 
values were used: 

- 1500 m/sec is the sound speed at the ocean surface and 
at the 2000 meter bottom, and 

- a minimum sound speed of 1475 m/sec occurs at the 1000 
meter depth. 

Three sets of simulations are presented using 5, 11, and 
17 sound-speed versus depth data pairs. The simulation 
results will be examined to quantify the number of ocean 
medium sound-speed samples required by each method to converge 
to a solution. Also, the agreement in results for the two 
methods and the CPU times will be evaluated. 

As the number of equally spaced data pairs varies, the 
medium and source conditions will be: 

- source depth of 1000 meters, 

- horizontal range of 3 kilometers, and 

five values of radial spatial frequency (FR) 
evaluated. 



27 



Tables 3.10, 3.11, and 3.12 show the results for the 5, 
11, and 17 data pair test cases for the piecewise linear 
method. In comparing the first two of these tables, large 
discrepancies are seen in the final depths for each of the 
five radial frequencies evaluated — errors as large as seven 
meters. Likewise, travel time values differ by hundredths of 
a second meaning that modulo 2 ir phase values show no 
correlation between test cases. Unfortunately, the variations 
between Tables 3.11 and 3.12 where more data points are used, 
show the piecewise linear method is not converging to a 
solution . 

Figures 3.7 through 3.9 correspond to Tables 3.10 through 
3.12. While numerical phase values show gross disagreement, 
the inability to detect differences in the graphical 
representations illustrates the very sensitive nature of the 
phase calculations. 

Tables 3.13 through 3.15 show the results when the ODE 
solver runs the simulations for 5, 11, and 17 data pairs. The 
5 data pair case of Table 3.13 again shows significant 
disagreement with the Table 3.14 results using 11 data pairs. 
Final depths show differences of as much as a meter while 
travel times vary by milliseconds. As with the piecewise 
linear method, the modulo 2n phase values have not yet 
converged to a solution. Variations of up to 2 radians are 
seen. 
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TABLE 3.10 

FIVE DATA PAIRS , PIECEWISE LINEAR SOLUTION 

RAY TRACING USING PIECEWISE LINEAR APPROXIMATION 
CASE 1 SOS 



NDATA * 5 NUMBER OF GRADIENTS » 4 



G( 1 ) » -0. 35355E-01 


1/SEC C( 


2 ) > -0. 14645E-01 


1/SEC 




G( 5) > 0.14645E-01 


1/SEC G ( 


4) ■ 0 . 55355E-01 


1/SEC 




ERAYS » F 




PRTALL • F 






YO = 1000.0 M 




DEPTH « 2000.0 M 


FREOC « 250.0 HZ 


CMAX 


FREQC/CMAX ■ 0 . 16667 


CYCLES/M 


RATIO » 0.9999 


FRMAX - 0.16665 CYCLES/M 


NFRS 


DLTFR » 0.035350 CYCLES/M 


RNGSTP > 1 .0 M 


HRZRNG » 3000.0 M 


YR > 



NOTE: FRM1N » DLTFR 



1500.0 M/SEC 
5 

<0.0 M 



NFR 


BETAOCDEG) 


DEPTH(M) 


TRVLT (SEC ) 


BETAY(DEG) 


THETAT(RAD) 


THETAR (RAD) 


THETAY(RAD) 


modulo twop: 

THETAY(RAD) 


1 


11.3<1 


144.74 


10.207186 


168.504 


0. 1603341E*05 


0 . 6282557E*03 


0. 1540515E*05 


5.0674943 


2 


23.160 


41.81 


5.103311 


156.449 


0.801 62 62E* 04 


0 . 12565 1 1E*04 


0 . 6 75975 1E*04 


5.3264822 


5 


36.153 


1064.46 


3.402439 


36. 180 


0 .5344538E *04 


0. 1884767E*04 


0.3459771E*04 


4*. 0191220 


*5 


51 .868 


682.47 


2.55474? 


127.902 


0. 401298 7E*04 


0.251 302 3E * 04 


0. 1499964E*04 


4 .5658515 


5 


79.49*5 


1510.35 


2.057830 


81.249 


0 . 3232433E*04 


0.31412 78E* 04 


0.91 15404E*02 


3. 1894409 



TOTAL CPU TIME » 0 MIN , 2*5.95 SEC 



TABLE 3.11 

ELEVEN DATA PAIRS , PIECEWISE LINEAR SOLUTION 



RAY TRACING USING PIECEWISE LINEAR APPROXIMATION 
CASE 1S11 



NDATA 


» 


11 




NUM8ER OF GRADIENTS * 


10 




G ( 1 ) 


> 


-0.38627E-01 


1/SEC 


G ( 2 ) - -0.34846E-01 


1/SEC 




G ( 3) 


• 


-0.27654E-01 


1/SEC 


G( 4 ) • -0. 17755E-01 


1/SEC 




G ( 5) 


* 


-0.61179E-02 


1/SEC 


G< 6) » 0.611 7 9E-02 


1/SEC 




G( 7) 


• 


0. 1 7755E-01 


1/SEC 


G( 8 ) = 0.27 654E-01 


1/SEC 




G( 9) 


* 


0. 34846E-01 


1/SEC 


GUO) » 0.38627E-01 


1/SEC 




ERAYS 


■ 


F 




PRTALL » F 






YO » 


1000.0 M 




DEPTH » 2000.0 M 


FREOC » 250.0 HZ 


CMAX « 1500.0 M/SEC 


FREOC/CMAX * 0.16667 


CYCLES/M RATIO * 0.9999 


FRMAX ■ 0.16665 CYCLES/M 


NFRS » 5 


DLTFR 


» 


0.033330 CYCLES/M 


RNGSTP » 1.0 M 


HRZRNG » 3000.0 M 


YR » 40.0 M 



NOTE: FRM1N » DLTFR 



NFR 


BETAOCDEG) 


DEPTH(M) 


TRVLT (SEC ) 


BETAY(DEG) THETAT(RAD) 


THETAR (RAD ) 


THETAY(RAD) 


MODULO TWOP 
THETAY(RAD) 


1 


1 1 . 34 1 


137.50 


10.216764 


168.506 0. 1604846E*05 


0 . 628255 7E *03 


0. 1542020E*05 


1.2628446 


2 


23.160 


37.97 


5.108073 


156.449 0 . 802 3 742E * 04 


0 . 12565 1 1E*04 


0. 6767231E*04 


0.2400057 


3 


36.153 


1067.39 


3.405626 


36.165 0 . 5 34 9545E* 04 


0. 1884767E*04 


0.34647 78E *04 


2.7430820 


4 


51 .868 


679.33 


2.557346 


127.965 0. 401707 OE* 04 


0.251 3023E * 04 


0 . 1504 048E 4 04 


2.3663691 


5 


79.494 


1520.40 


2.060168 


81.354 0 . 32 36 1 04E * 04 


0.3141278E*04 


0 . 9482580E *02 


0.5780165 


TOTAL CPU TIME 


» 0 MIN 


. 25.36 SEC 
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TABLE 3.12 

SEVENTEEN DATA PAIRS , PIECEWISE LINEAR SOLUTION 



RAY TRACING USING PIECEWISE LINEAR APPROX 1 MAT 1 ON 
CASE IS 1 7 



NDATA 


a 


17 




NUMBER 


OF GRADIENTS a 


16 




G ( 1) 


a 


-0 . 39018E-0 1 


1/SEC 


GC 2) 


» -0.3751 9E-0 1 


1/SEC 




G ( 3) 


a 


-0 . 345 77E-0 1 


1/SEC 


G ( 4 ) 


« -0.3030 7E- 0 1 


1/SEC 




G ( 5 ) 


a 


-0 .24873E-01 


1/SEC 


G ( 6) 


a -0. 18482E-01 


1/SEC 




G ( 7) 


a 


-0.1138 1E-0 1 


1/SEC 


G ( 8) 


■ -0.38429E-02 


1/SEC 




G ( 9) 


■ 


0 . 38 429E- 02 


1/SEC 


G( 10) 


a 0.11381E-01 


i/SEC 




G( 11 ) 


a 


0. 18482E-01 


1/SEC 


G( 1 2 ) 


a 0. 24873 E- 01 


1/SEC 




G ( 1 3 ) 


a 


0 .30307E-01 


1/SEC 


G< ]< ) 


■ 0 . 345 7 7E-0 1 


1/SEC 




G C 1 5 ) 


a 


0. 37519E-01 


1/SEC 


G ( 1 6 ) 


» 0. 39018 E- 01 


1/SEC 




ERAYS 


a 


F 




PRTALL » F 






YO « lOOO.O M 




DEPTH > 2000.0 M 


FREOC a 250.0 HZ 


CMAX a 1500.0 M/SEC 


FREOC/CMAX « 0.16667 


CYCLES/M RATIO » 0.9999 


FRMAX » 0.16665 CYCLES/M 


NFRS ■ 5 


DLTFR 


a 


0.033330 CYCLES/M 


RNGSTP a 1.0 M 


HRZRNG « 3000.0 M 


YR ■ 40.0 M 


NOTE: 


FRM1N » DLTFR 













NFR 


BETA0 ( DEG ) 


DEPTH(M) 


TRVLT ( SEC ) 


BETAY ( DEG ) 


THETAT (RAD ) 


THETAR(RAD) 


THETAY(RAD) 


MODULO TWOPI 
THETAY(RAD) 


1 


11.341 


136.67 


10.217866 


168.506 


0. 1605019E*05 


0.6282557E*03 


0 . 15421 93E*05 


2.9937577 


2 


23.160 


37.53 


5.108620 


156.449 


0.8024602E*04 


0.125651 1E*04 


0. 676809 1E*04 


1.1002238 


3 


36.153 


1067.73 


3.405996 


36.160 


0.5350126E*04 


0 . 1884767E*04 


0 . 346535 9E*04 


3.3242727 


4 


51.868 


678.98 


2.557644 


127.973 


0.401 7538E*04 


0.251 3023E*04 


0.1504515E*04 


2.8338755 


5 


79.494 


1521.62 


2.060452 


81.322 


0 . 3236550E*04 


0 . 3141278E*04 


0. 9527150E 4 02 


1.0237247 


TOTAL CPU TIME ■ 0 MIN 


. 25.94 SEC 
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f 



Figure 3.7 Ray trace corresponding to Table 3.10 
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Figure 3.8 



Ray trace corresponding to Table 3.11 
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ray trace:: piecewise: linear approximation 




Figure 3.9 Ray trace corresponding to Table 3.12 



33 



CASE: 1S17 

YO: 1000.0 M HRZRNG: 3.0 KM 

RAYS: 5 



TABLE 3.13 

FIVE DATA PAIRS, ODE SOLVER 

RAY TRACING USING AK1MA CU8IC SPLINE « ODE SOLVER 
CASE 2S05 

NDATA - S 

ERAYS « F PR TALL • F 



Y0 • 1000.0 M 


DEPTH • 2000.0 H 


FREOC « 


250.0 H2 


CMAX • 


1500.0 M/SEC 


FREOC/CMAX . 0.16667 CYCLES/H 


RATIO ■ 0.9999 


FRHAX « 0. 


1 6665 CYCLES/H 


NFRS • 


5 


DLTFR • 0.035550 CYCLES/M 


RMGSTP - 1.0 M 


HR2RNG ■ 


5000.0 M 


YR - 


<0.0 M 



NOTE J FRMIN ■ DLTFR 



NFR 


BETAO(DEG) 


DEPTM(M) 


TRVLT (SEC ) 


BETAYC DEG 1 


THETATCRAD 1 


THET AR (RAD 1 


THET AY (RADI 


MODULO TMOP] 
THETAY(RAD) 


1 


1 1 . 5< 1 


1 5< .85 


10.220257 


168.509 


0 . 1 <0539<E<05 


0 . 6282557E<03 


0.1S<2S69E<05 


0.4666177 


2 


25. U0 


36.56 


5 . 1 0 °8 2 3 


156. <52 


0 . 8026<91E<0< 


0 . 125651 1E<0< 


0 . 6 7 6 990 0E<04 


2.9891864 


3 


56.153 


1068. <9 


3 .<06808 


36.157 


0. 53SI<02E<0< 


0. 18847<7E<0< 


0 . 3<66635E<0< 


4.5997015 


< 


51.868 


678.30 


2 .558 194 


127.973 


0. <018402E<0< 


0.251 3023E*0< 


0 . 1505 3 7 9E <04 


5.6978733 


5 


7 9 . < 9< 


1521.82 


2 . 0 6 0< 98 


81.306 


0. 323662 3E *04 


0. 31<1278E<0< 


0 . 9534<<8E<02 


1.0968988 



TABLE 3 . 14 

ELEVEN DATA PAIRS, ODE SOLVER 

RAY TRAC1N0 USING AK IMA CUBIC SPLINE t ODE SOLVER 
CASE 2 SI 1 

NDATA • 11 



ERAYS • F 


PR TALL « F 






Y0 • 1000.0 M 


DEPTH - 2000.0 M 


FREOC • 250.0 H2 


CMAX • 1500.0 M/SEC 


FREOC/CMA* - 0.16667 CYCLES/M 


RATIO » 0.9999 


FRMAX - 0.16665 CYCLES/M 


NFRS - 5 


DLTFR • 0.033330 CYCLES/M 

NOTE 1 FRMIN • DLTFR 


RNGSTP - 1.0 M 


HRZRNG ■ 3000.0 H 


YR ■ 40.0 M 



NFR 


BETA0(DEG 1 


DEPTH CM 1 


TRVLT ( SEC ) 


BETAY(DEG 1 


THETAT (RAD) 


TRETAR (RAD ) 


THETAY(RAD) 


MODULO TMOP] 
THETAY(RAD) 


1 


11.341 


136.09 


10.218623 


168.506 


0 . 1 6051 38E<05 


0 . 628255 7E<03 


0 . 1542312E<05 


4.1828982 


2 


23.160 


37.23 


5.108997 


156. «9 


0. 802519<E<0< 


0 . 125651 1E<0< 


0. <768683E<0< 


1 .6920625 


3 


36.153 


1067.97 


3.406254 


36.157 


0 . 5 3505 3 IE *04 


0. 18847<7E<0< 


0. 346576<E<0< 


3.7288771 


4 


51.868 


6 78 .73 


2.557848 


127.978 


0 .4017858E*0< 


0 . 25 1 3023E* 04 


0. 150<835E*04 


3.1540595 


5 


79.494 


1522.39 


2.060633 


81.312 


0 . 32 3 6834E<04 


0. 31<1278E<0< 


0 . 9555559E*02 


1.3078080 



TOTAL CPU TIKE « <8 MIN , 14.53 SEC 



TABLE 3.15 

SEVENTEEN DATA PAIRS, ODE SOLVER 

RAY TRACING USING AKlHA CU8IC SPLINE 8 ODE SOLVER 
CASE 2S1 7 

NDATA - 17 

ERAYS • F PR TALL - F 

YO • 1000.0 M DEPTH - 2000.0 M FREOC « 250.0 M2 CMAX • 1500.0 M/SEC 

FREOC/CHAX • 0.16667 CYCLES/M RATIO ■ 0.9999 FRHAX - 0.16665 CYCLES/H NFRS « 5 



DLTFR 

NOTES 


- 0.033330 CYCLES/H RNGSTP ■ I. 

FRMIN - DLTFR 


0 M HR2RNG 


- 3000.0 M 


YR ■ 40.0 


i H 


NFR 8ETA0 ( DEG 1 


DEP TH(M 1 


TRVLT (SEC ) 


8ETAY(DEG ) 


THETAT (RAD ) 


THETAR (RADI 


THETAY(RAD) 


MODULO TWOP1 
TMETAY(RAD 1 


1 


11.341 


136.13 


10.218580 


160.506 


0. 1605131E<05 


0 . 6282557E<03 


0.1542305E*05 


4.1156874 


2 


23.160 


37.25 


5.108975 


156.449 


0 .8025 159E<04 


0. 125<511E<0< 


0. 67686<8E<0< 


1.6570617 


3 


36.153 


1067.96 


3.406239 


36.157 


0 .5550507E*04 


0. 188<767E<0< 


0. 3<657<0E<04 


3.7048568 


< 


51 .868 


6 78 .74 


2.557036 


127.978 


0.4017839E<04 


0.2513023E*04 


0. 150<817E<0< 


3.1353514 


5 


79.494 


1522.36 


2.060625 


01.317 


0. 323<822E<0< 


0. 31<1278E<0< 


0. 955<37<E<02 


1 .2959589 



TOTAL CPU TIME ■ <7 MIN , <5.20 SEC 
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A comparison of Tables 3.14 and 3.15 shows that phase 
calculation results converged for each of the five FR values. 
The increase in data pairs to 17 from 11 has changed the 
final depth values by only a few hundredths of a meter. 
Likewise, the very sensitive modulo 2Tt phase values show a 
change of only hundredths of a radian or less. The ray trace 
plots for Tables 3.13 through 3.15 are provided in Figures 
3.10 through 3.12. 

The CPU times for these three tables show the ODE solver 
to be costly, but insensitive to increases in the number of 
data pairs. That is, approximately the same amount of CPU 
time is required regardless of the number of data pairs used. 
The CPU time for the piecewise linear method also proved 
fairly insensitive to the number of data pairs used as seen in 
Tables 3.10 through 3.12. 

The significant findings of this section are: 

- the ODE solver can perform accurate phase 
calculations with only 17 data pairs sampled from a 
2000 meter deep ocean having a curved, sinusoidal 
sound-speed profile, 

- the piecewise linear approach does not converge to a 
solution using the 17 data pairs, and 

- the differential equation solver is very costly to use 

in terms of CPU time versus the piecewise linear approach. 
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Figure 3.10 



Ray trace 



corresponding 



to Table 3.13 
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RRY TRACE:: rkihr cubic spline: l ODE SOLVER 
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Figure 3.11 



Ray trace corresponding to Table 3.14 
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HRZRNG: 3.0 KM 



RflY TRACE: RKIMR CUBIC SPLINE & ODE SOLVER 




Figure 3.12 Ray trace corresponding to Table 3.15 
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CASE: 2S17 

YO: 1000.0 M HRZRNG: 3.0 KM 

RAYS: 5 




E. SPECIAL INVESTIGATIONS 



Each of the two phase calculation methods has shown one 
characteristic weakness. The piecewise linear method did not 
converge to a solution in the sinusoidal sound-speed profile 
case when 17 data pairs or samples were provided. In 
contrast, the ODE solver does converge to a solution, but is 
costly in terms of CPU time. This section will attempt to 
tailor the settings for each method to compensate for its 
weaknesses . 

The simulation results in Tables 3.16 through 3.18 were 
produced by the piecewise linear method using 29, 55, and 65 
data pairs, respectively. All medium and sound-speed profile 
characteristics remain the same as in the previous section. 
As the number of data pairs increases, the 0 y values are 
converging to the Table 3.15 solutions obtained with the ODE 
solver. Unfortunately, the phase values continue to show 
unacceptable errors of tenths and hundredths of a radian for 
the best case of 65 data pairs. 

Table 3.19 shows the results for the ODE solver when the 
range step (RNGSTP) parameter is increased to five meters. As 
explained earlier, the range step is the differential 
equation system independent variable. All calculations prior 
to Table 3.19 used a conservative range step of one meter for 
both methods. Table 3.19 compares favorably with Table 3.15 
achieving a balance between CPU time and accuracy. While the 
degradation in accuracy is only microseconds for travel time 



39 



TABLE 3.16 

TWENTY-NINE DATA PAIRS , PIECEWISE LINEAR SOLUTION 



RAY TRACING USING PIECEWISE LINEAR APPROXIMATION 
CASE 1S29 



NDATA 


■ 


29 




NUMBER 


OF GRADIENTS ■ 


» 28 


G ( 1) 


■ 


-0.39188E-01 


1/SEC 


G ( 2) 




-0.38 6 95E- 0 1 


1/SEC 


G ( 3) 


■ 


-0.37715E-01 


1/SEC 


G ( 4) 




-0 . 36262E- 0 1 


1/SEC 


G ( 5) 


■ 


-0.34352E-01 


1/SEC 


G ( 6) 




-0.3201 OE-O 1 


1/SEC 


G( 7) 


■ 


-0.29266E-01 


1/SEC 


G ( 8) 




-0.26 154E- 0 1 


1/SEC 


G ( 9) 


• 


-0.22712E-01 


1/SEC 


G ( 1 0 ) 




-0. 18986E-01 


1/SEC 


GUI) 


• 


-0. 15020E-01 


1/SEC 


G ( 1 2 ) 




-0 . 1 0866E-01 


1/SEC 


G ( 1 3 ) 


• 


-0 . 65 745E-02 


1/SEC 


G ( 1 4 ) 




-0 . 22007E-02 


1/SEC 


G( 15 ) 


• 


0.22007E-02 


1/SEC 


G( 1 6 ) 




0. 65745 E- 02 


1/SEC 


G ( 1 7 ) 


• 


0.10866E-01 


1/SEC 


G( 18 ) 




0. 15020E-01 


1/SEC 


G( 1 9 ) 


• 


0. 18986E-01 


1/SEC 


G ( 2 0 ) 




0.2271 2E- 0 1 


1/SEC 


G ( 2 1 ) 


■ 


0. 26 154 E- 01 


1/SEC 


G ( 22 ) 




0 . 29266E-01 


1/SEC 


G (2 3 ) 


• 


0.3201 0E-01 


1/SEC 


G ( 24 ) 




0. 34352E-01 


1/SEC 


0(25) 


■ 


0. 36262 E- 01 


1/SEC 


G (2 6 ) 




0.377 15E-01 


1/SEC 


G(27) 


■ 


0 • 38495E-01 


1/SEC 


G ( 28 ) 




0.39 188E-0 1 


1/SEC 


ERAYS 


• 


F 




PRTALL ■ F 




Y0 • 


1000.0 M 




DEPTH ■ 2000.0 M 


FRI 



FREOC/CMAX ■ 0.16*47 CYCLES/M RATIO ■ 0.9999 FRMAX - 
DLTFR • 0,033330 CYCLES/H RNGSTP ■ 1,0 M HR2RN0 
NOTEi FRM1N • DLTFR 



250.0 HZ 

0.16 665 CYCl.ES/M 
» 3000.0 M 



CMAX • 1500.0 M/SEC 

NFRS ■ 5 

YR ■ 40.0 M 



NFR 


BETAO(DEG) 


DEPTH(M) 


TRVLT(SEC) 


BETAY(DEG) 


THETATCRAD) 


THETAR (RAD) 


THETAY(RAD) 


MODULO twop: 
THETAYCRAD) 


1 


11.341 


134.31 


10.218342 


168.506 


0 . 1605093E 4 05 


0 . 6282557E*03 


0.1542248E»05 


3.7413092 


2 


23.160 


37.34 


5.108857 


154.449 


0 .8024 97 3E*04 


0. 12565 1 1E*04 


0 . 6 748462E*04 


1.4710844 


3 


36.153 


1067 .88 


3.406158 


36.157 


0 .5350380E*04 


0 . 1884 7 6 7E ♦ 04 


0 • 34656 1 3E*04 


3.5775983 


4 


51 .868 


678.82 


2.557772 


127.976 


0 . 4 0 1 7 7 3 9E*04 


0 .25 1 3023E*04 


0. 1504716E»04 


3.0348743 


5 


79.494 


1522.11 


2.060567 


81.320 


0 . 32 34 730E*04 


0.314 12 78E*04 


0 . 9545 1 97E*02 


1.2041934 



TOTAL CPU TIME • 0 MIN , 26.46 SEC 
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TABLE 3.17 

FIFTY-FIVE DATA PAIRS , PIECEWISE LINEAR SOLUTION 



RAY TRACING USING PIECEWISE LINEAR APPROXIMATION 
CASE 1S55 



NDATA 




55 




NUMBER 


OF GRADIENTS » 


54 


g ( n 




-0.39248E-01 


i/SEC 


G( 2) 




-0.391 15E-0 1 


l/SEC 


G ( 3) 




-0 . 3S850E-0 1 


i/SEC 


G ( 4) 




-0.3845 3E-0 1 


l/SEC 


G< 5) 




-0.37 926E-0 1 


1/SEC 


G ( 6) 




-0 . 37271E-01 


l/SEC 


G( 7) 




-0. 36490E-01 


i/SEC 


G( 8) 




-0.355S6E-01 


l/SEC 


Gt 9) 




-0.34561E-01 


1/SEC 


G ( 1 0 ) 




-0.334 1 9E-0 1 


l/SEC 


G(1I) 




- 0 . 32 1 6 3E- 0 1 


i/SEC 


G ( 12 ) 




-0 .30800E-01 


l/SEC 


G( 13 ) 




-0.29331E-01 


i/SEC 


G ( 1 4 ) 




-0.27764E-01 


l/SEC 


G ( 15 ) 




-0.2610 3E-0 1 


1/SEC 


G ( 1 6 ) 




-0 . 2435 3E-0 1 


l/SEC 


G( 1 7 ) 




-0.22521E-01 


1/SEC 


G ( 18 ) 




-0.20413E-01 


l/SEC 


G(19) 




-0.18 6 35E-0 1 


1/SEC 


G(20 ) 




-0. 16594E-01 


l/SEC 


GC21) 




-0.14497E-01 


1/SEC 


G (22 ) 




-0. 12350E-01 


l/SEC 


G ( 2 5 ) 




-0.101 62E-0 1 


1/SEC 


G (24 ) 




-0. 79399E-02 


l/SEC 


G (25 ) 




-0.56907E-02 


1/SEC 


G (26 ) 




-0 . 3422 IE-02 


l/SEC 


0(27) 




-0 . 11420E-02 


1/SEC 


G ( 28 ) 




0.11420E-02 


l/SEC 


G(29) 




0. 34221 E- 02 


i/SEC 


G ( 30 ) 




0 . 56907E-02 


l/SEC 


G ( 3 1 ) 




0 . 79399E-02 


1/SEC 


G ( 32 ) 




0.101 62E-0 1 


l/SEC 


0(33) 




0 . 12350E-01 


i/SEC 


G ( 34 ) 




0 . 14497E-0 1 


l/SEC 


G( 35 ) 




0 . 1 6594E-0 1 


1/SEC 


G ( 36 ) 




0. 19635E-01 


l/SEC 


G( 37 ) 




6 . 2061 3E-01 


1/SEC 


G( 38 ) 




0 . 2252 1E-0 1 


l/SEC 


0(39) 




0.24353E-01 


1/SEC 


G ( 4 0 ) 




0.261 0 3E-0 1 


l/SEC 


G(41) 




0.27764E-01 


1/SEC 


G(42) 




0. 29331 E- 01 


l/SEC 


G(43) 




0 . 30800E-0 1 


1/SEC 


G (44 ) 




0 . 32 163E-01 


l/SEC 


G(45) 




0.3341 9E-0 1 


l/SEC 


0(46) 




0. 34561E-01 


l/SEC 


G(47) 




0 * 35586E-0 1 


1/SEC 


G ( 48 ) 




0.364 90E-01 


l/SEC 


G(49) 




0. 37271 E- 01 


l/SEC 


G(50) 




0 . 37926E-0 1 


l/SEC 


G(51 ) 




0 . 38453E-0 1 


l/SEC 


G(52 ) 




0 . 38850E-0 1 


l/SEC 


G(53) 




0.391 15E-0 1 


l/SEC 


G ( 54 ) 




0 . 39248E-01 


l/SEC 



ERAYS • F 


PRTALL ■ P 








Y0 • 1000.0 M 


DEPTH • 2000.0 M 


FREOC • 


250.0 HZ 


CMAX • 1500.0 M/SEC 


FREOC/CMAX • 0.16667 CYCLES/M 


RATIO • 0.9999 


FRMAX ■ 0. 


,16665 CYCLES/M 


NFRS • 5 


DLTFR • 0.033330 CYCLES/M 

NOTE: FRM1N ■ DLTFR 


RNGSTP • 1.0 M 


HRZRNG • 


3000. 0 M 


YR • 40.0 M 



NFR 


BETAO(DEG) 


DEPTH(M) 


TRVLT (SEC) 


BETAY(DEG) 


THETAT (RAD ) 


THETAR(RAD) 


THETAY(RAD) 


MODULO TWOPl 
THETAY(RAD) 


1 


11.341 


136.18 


10.218510 


168.506 


0 . 1 605 120E*05 


0 . 6282557E*03 


0.1542294E*05 


4.0060532 


2 


23.160 


37.28 


5.108940 


156.449 


0 .80251 04E 4 04 


0 . 125651 1E+04 


0 . 6768593E*04 


1,6024783 


3 


36.153 


1067.94 


3.406215 


36.157 


0.5350470E*04 


0.1884767E404 


0 . 3465702E*04 


. 3.6673828 


4 


51.868 


678.77 


2.557817 


127.977 


0 . 40 1 78 1 0E* 04 


0.2513023E*04 


0. 1504787E + 04 


3.1058334 


5 


79.494 


1522.29 


2.060608 


•1.317 


0 . 32367 95E*04 


0.31412 78E 4 04 


0 .955 1 684E 4 02 


1.2690645 



TOTAL CPU TIME • 0 MIN , 26.61 SEC 
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TABLE 3.18 

SIXTY-FIVE DATA PAIRS, PIECEWISE LINEAR SOLUTION 



RAY TRACING USING PIECEWISE LINEAR APPROXIMATION 
CASE 1 S65 



NDATA 




65 




NUMBER OF GRADIENTS ■ 


6 A 


G ( 1 ) 




-0. 3925AE-01 


1/SEC 


3 ( 2) 




-0. 39160E-01 


1/SEC 


G( 3) 




-0. 38971E-01 


1/SEC 


G ( A) 




-0. 38688E-01 


1 /SEC 


G ( S) 




-0 .383 12E-01 


1/SEC 


G ( 6) 




-0 . 378AAE-01 


1/SEC 


G( 7) 




-0.3728AE-01 


1/SEC 


G( 8) 




-0. 36635E-01 


1/SEC 


G ( 9 ) 




-0 . 35897E-01 


1/SEC 


GUO) 




-0 . 35073E-0 1 


i/SEC 


GUI) 




-0.3A165E-01 


1/SEC 


G( 1 2 ) 




-0. 3317AE-01 


1/SEC 


G ( 1 3 ) 




-0.321 0 3E- 0 1 


1/SEC 


G ( 1A ) 




-0 . 30955E-01 


1/SEC 


G C 1 5 ) 




-0.29733E-01 


1/SEC 


G( 16) 




-0 . 28 A 38E- 0 1 


1/SEC 


G ( 1 7 ) 




-0.2 7075E- 01 


1 /SEC 


G( 18 ) 




-0.256A7E-01 


1/SEC 


G ( 1 9 ) 




- 0 . 2 A 1 58E- 0 1 


1/SEC 


G ( 2 0 ) 




-0.2261 0E- 0 1 


1/SEC 


G ( 2 1 ) 




-0.21007E-01 


1/SEC 


G(22 ) 




-0. 1935AE-01 


1/SEC 


G ( 23 ) 




-0. 1 765AE-01 


1/SEC 


G(2A ) 




-0. 15912E-01 


1/SEC 


G ( 25 ) 




-0. 1A 132 E- 01 


1/SEC 


GC 26 ) 




-0.1231 7E-0 1 


l/SEC 


G ( 2 7 ) 




-0. 10A73E-01 


1/SEC 


G C 28 ) 




-0 . 8 6032E-02 


1/SEC 


G ( 2 9 ) 




-0.671 30E-02 


1/SEC 


G( 30 ) 




-O.A8066E-02 


l/SEC 


G ( 3 1 ) 




-0.28886E-02 


1/SEC 


G( 32 ) 




-0 . 9636AE-03 


l/SEC 


G ( 3 3 ) 




0.9636AE-03 


1/SEC 


GC3A ) 




0 . 2888 6E-02 


l/SEC 


G ( 35 ) 




0.A8O66E-02 


1/SEC 


G ( 3 6 ) 




0.671 30E- 02 


l/SEC 


G ( 3 7 ) 




0 . 8 6032E-02 


1/SEC 


G(!8) 




0. 10A73E-01 


l/SEC 


G ( 3 9 ) 




0. 12S17E-01 


1/SEC 


G ( A 0 ) 




0. 1A132E-01 


l/SEC 


cun 




0. 15912E-01 


1/SEC 


G(A2) 




0. 1765AE-01 


l/SEC 


G ( A 3 ) 




0. 1935AE-01 


1/SEC 


G ( A A ) 




0.21 00 7E- 0 1 


l/SEC 


G(A5 ) 




0.2261 OE-O 1 


1/SEC 


G ( A 6 ) 




0.2A158E-01 


l/SEC 


q«7) 




0.256A7E-01 


1/SEC 


G(A8) 




0.270 75E-0 1 


l/SEC 


G C A 9 ) 




0 . 28 A 38E-0 1 


1/SEC 


G(50 ) 




0.297 3 3E-0 1 


l/SEC 


G ( 5 1 ) 




0.30955E-01 


1/SEC 


G( 52 ) 




0.32103E-01 


l/SEC 


G (5 3 ) 




0.3317AE-O1 


1/SEC 


G(5A ) 




0. 3A165E-01 


l/SEC 


G( 55 ) 




0. 35073E-01 


1/SEC 


G( 5 6 ) 




0.358 9 7E-0 1 


l/SEC 


G( 5 7 ) 




0 . 36635E-01 


1/SEC 


G ( 58 ) 




0.3728AE-01 


l/SEC 


G(5 9 ) 




0 . 3 78 A AE-0 1 


1/SEC 


G ( 6 0 ) 




0. 38312E-01 


l/SEC 


G( 61 ) 




0. 38 688E-01 


1/SEC 


G ( 62 ) 




0 . 38971E-01 


l/SEC 


G ( 6 3 ) 




0. 39160E-01 


1/SEC 


G(6A ) 




0.3925AE-01 


l/SEC 



ERAYS « F 


PR TALL = F 






YO » 1000. 0 M 


DEPTH ■ 2000.0 M 


FREOC ■ 250.0 HZ 


CMAX - 1500.0 H/SEC 


FREOC/CMAX - 0.16667 CYCLES/M 


RATIO « 0.9999 


FRMAX « 0.16665 CYCLES/H 


NFRS « 5 


DLTFR • 0.033330 CYCLES/M 

NOTE: FRM1N - DLTFR 


RNGSTP » 1 .0 M 


HRZRNG • 3000. 0 M 


YR ■ A 0 . 0 H 



NFR 


BETAO(DEG) 


DEPTH(M) 


TRVLT ( SEC ) 


BETAY(DEG) 


THETAT ( RAD ) 


THETAR (RAD ) 


THETAY(RAD) 


MODULO TWOP 
THETAY (RAD ) 


1 


11 . 3A1 


136.17 


10.218528 


168.506 


0. 1 605 123E*05 


0 628 25 5 7E * 0 3 


0 . 15A2297E* 05 


A.03A0889 


2 


23.160 


37.27 


5 . 1 08 9 A 9 


156. AA9 


0 .8025 1 18E* OA 


0.125651 1E*0A 


0 . 6 768 607E 4 OA 


1 .6 163983 


3 


36.153 


1067 .9A 


3 . A 0622 1 


36.157 


0.5350A79E«OA 


0.188A767E*OA 


0 . 3A65 7 12E 4 OA 


3.6769515 


A 


51 .868 


678 . 76 


2.557822 


127 .978 


0.A017818E*0A 


0 .25 1 3023E*0A 


0. 150A795E 4 0A 


3.113A251 


5 


79.A9A 


1522.31 


2.060612 


81.317 


0. 3236802E‘0A 


0.31A1278E*0A 


0 . 9552 3 6 OE* 02 


1 .2758237 


TOTAL CPU TIME 


» 0 MIN 


, 26. 8A SEC 
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TABLE 3.19 

INCREASED RANGE STEP, ODE SOLVER SOLUTION 



RAY TRACING USING AKIMA CUBIC SPLINE i ODE SOLVER 
CASE 2S17 

NDATA - 17 



ERAYS » F 
YO ■ 1000.0 M 

FREOC/CMAX • 0.16667 CYCLES/M 



PRTALL ■ F 
DEPTH - 2000.0 H 
RATIO * 0.9999 



FREOC 

FRMAX 



250.0 HZ 

0. 16665 CYCLES/M 



CMAX 

NFRS 



1500.0 M/SEC 
5 



DLTFR 


» 0.033330 CYCLES/M RNGSTP - 5. 


0 M HRZRNG 


• 3000.0 M 


YR * 40.0 


M 


NOTE: 


FRM1N • DLTFR 














NFR BETA0CDEG) 


DEPTHCM) 


TRVLTCSEC) 


BETAY ( DEG ) 


THETAT (RAD ) 


THETAR(RAD) 


THETAY(RAD) 


MODULO TWOPJ 


















THETAYCRAD) 


1 


11.341 


136.09 


10.218584 


168.506 


0. 1 6 OS 1 3 1 E ♦ OS 


0.6282557E*03 


0. 1542306E*05 


4.1219846 


2 


23.160 


37.25 


5 . 108 9 75 


156.449 


0 . 8 025 1 5 9E* 04 


0. 125451 1E*04 


0.47 68 648E * 04 


1.4570457 


3 


36.153 


1067 . 96 


3.406239 


36.157 


0 .5350507E*04 


0 . 18847 67E*04 


0.346S740E*04 


3.7047770 


4 


51 .868 


678.74 


2.557836 


127. 978 


0 . 40 1 78 39E ♦ 04 


0.251 S023E+04 


0. 1504817E*04 


3.1353500 


5 


79.494 


1522.36 


2.060625 


81.317 


0.3236822E*04 


0 . 3 1 4 1 2 78E *04 


0 . 9554 379E*02 


1 .2960153 


TOTAL 


CPU TIME 


• 10 MIN 


. 3.53 SEC 













and hundreds or thousandths of a radian for modulo 27 r 0 Y , the 
reduction in CPU time is following a nearly linear relation 
to the increase in range step, that is, CPU time was reduced 
approximately by one-fifth. 



F. EIGENRAYS 

Tables 3.20 and 3.21 show the eigenrays found by the two 
ray acoustics methods. Figures 3.13 and 3.14 correspond to 
these two tables. The figures show the sound-speed profile to 
be that of a SOFAR channel. This channel profile was chosen 
because it produces interesting ray patterns, and it commonly 
occurs in nature. The minimum sound-speed occurs at 1000 
meters depth and a local minima occurs at the ocean surface. 
The tables show that each method found the same eigenrays at 
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TABLE 3.20 

EIGENRAYS FOUND USING THE 
PIECEWISE LINEAR APPROXIMATION METHOD 



RAV TRACING USING PIECEWISE LINEAR APPROXIMATION 
CASE 1CM4 



NDATA = 4 

G( 1) = 0. 1 6 OOOE-O I l/SEC 

G ( 3) = 0 . 1 7000E-0 1 l/SEC 



NUH8ER OF GRADIENTS = 3 

G ( 2 1 * -0 . 1888 9E- 0 1 l/SEC 



ERAYS = T 
YO a 100.0 H 
RNGSTP » 5.0 M 

YERROR * 5.0 H 

RAY BETAO(DEG) DEPTH(M) 



PR TALL = F 
DEPTH * 2000.0 H 
ANGSTP = 1.0 DEG 



TRVLT(SEC) 8ET AY( DEG ) 



FREOC = 250.0 HZ 

HRZRNG « 200.0 H 

EIGENRAYS 

THETAT(RAD) THETAR(RAD) 



YR = 250.0 



THETAY ( R AD ) 



1 53.000 251.11 

2 54.000 245.69 

3 150.000 246.73 



0.167269 52.856 
0.165112 53.856 
0.267040 29.939 



0 . 262 74 62E *0 3 
0 . 25 935 6 ?E 4 0 3 
0 .4 194649E*03 



0. 1 672491E*03 
0. 1694232E*03 
0. 10470?3E*03 



0 .95497 10E 4 02 
0.899337 IE *02 
0.3147556E*03 



TOTAL CPU TIME * 0 H1N , 22.47 SEC 



TABLE 3.21 

EIGENRAYS FOUND USING THE 
ODE SOLVER METHOD 



RAY TRACING USING AK IMA CUBIC SPLINE * ODE SOLVER 
CASE 2CH4 

NDATA « 4 

YR = 250.1 

TMETAY(RAD) 

0 . 95 1 7250E* 02 
0 . 8 °620 32E 4 02 
0 . 3 14447OE*03 

TOTAL CPU T 1 HE = 24 H1N , 12.94 SEC 



ERAYS 



PRTALL = F 



DEPTH = 2000.0 H 
ANGSTP * 1.0 DEG 



FREOC * 
HRZRNG « 



YO * 100.0 H 

RNGSTP * 5.0 H 

YERROR s 5.0 M 

EIGENRAYS 

RAY BET AO ( DEG ) DEPTH(H) TRVLT ( SEC ) BETAYt DEG I THE T AT ( R AD ) 



250.0 MZ 
200.0 H 



THETAR (RAD ) 



1 53.000 250.85 0.167063 52.930 

2 54.000 245.44 0.164912 53.932 

3 150.000 246.56 0.266843 29.971 



0. 262421 6E*03 
0 .25 <j 0435E»03 
0. 419156 3E*03 



0 . 1 67249 IE*03 
0. 1694232E*03 
O.1047O93E*03 



H 



HODULO TWOP1 
THETAY (RAD) 

1 .2493187 
1.9691119 
0.5963779 



M 



HODULO TWOPI 
THETAY(RAD) 
0.924723? 
1.6557224 
0.2877732 
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RRY TRACE: P'ECEHISE LINEPR APPROX I MAT I OH 




Ray trace plot corresponding 
4 5 



Figure 3.13 



to Table 3.20 



MOO 1500 1600 0.00 0.02 0.01 0.06 0.03 0.10 0.12 0.M 0.16 

Cm/SJ RRNGE(KM) 

CASE: 1CH4 

YO: 100.0 M HRZRNG: 0.2 Ktt YR : 250.0 M 

EIGENRflYS: 3 YERROR: 5.0 PI ___ 




SSP RRY TRACE: flKIMQ CUBIC SPLINE i OOE SOLVER 




Figure 3.14 Ray trace plot corresponding to Table 3.21 
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launch angles of 53, 54 and 150 degrees. The tables highlight 
this problem's parameters to be: 

- a source depth y 0 of 100 meters; 

- a receiver depth y r of 250 meters; 

- and a horizontal range (HRZRNG) to the receiver of 200 
meters . 

As expected, the ODE solver runs into range limitations 
imposed by CPU time requirements. The 200 meter example 
presented consumed over 24 minutes of CPU time using the ODE 
solver versus the 21 seconds used by the piecewise linear 
method. 

Tables 3.22 and 3.23 (along with their corresponding 
Figures 3.15 and 3.16) show solutions to long-range eigenray 
search problems. These two examples continue to use a source 
depth y 0 of 100 meters and a receiver depth y r of 250 meters. 
Table 3.22 shows a horizontal range (HRZRNG) of 10 kilometers 
with a depth error y error at the receiver of 2 meters. The 
angle step size (angstp) between each ray is 0.1°. Table 3.22 
also shows that three eigenrays were found with launch angles 
of 51.2°, 90.9° and 129.5°. The CPU time for this 10km run is 
shown to be 8.5 minutes. 
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LONG RANGE 



LINEAR SOLUTION 



TABLE 3.22 

EIGENRAYS, PIECEWISE 



RAY TRACING USING PIECEWISE LINEAR APPROXIMATION 
CASE 1CHA 



NDATA = A 

G( 1) = 0. 1 60 OOE-O 1 1/SEC 

G( 3) * 0 . 1 7000E-01 l/SEC 



NUMPER OF GRADIENTS = 3 

G ( 2 ) - -0. 1888 fl E-01 l/SEC 



ERAYS = T 
Y0 = 100.0 M 

RNGSTP * 5.0 M 

YERRDR = 2.0 M 



PRTALL = F 
DEPTH = 2000.0 M 
ANGSTP = 0.1 DEG 



RAY BETA0 (DEG ) DEPTH (M) TRVLT (SEC ) BETAY(DEG) 



1 51.200 

2 90.900 

3 129.500 



2A8.91 
250.77 
2 A 9 . 2 1 



8.6A5958 

6.67A996 

8.730378 



51.067 

86.356 

50.370 



FREOC = 
HRZRNG 

EIGENRAYS 

THETAT(RAD) 

0 . 1 358 10AE*05 
0. 10A8506E*05 
0 . 1 3 7 1 3 65E* 05 



250.0 HZ 
10000.0 M 



THE TAR (RAD ) 

0 . 8 1 60 3 92E* 0A 
0 . 1 0A696AE*05 
0 .8079626E*0A 



YR = 250.0 M 



THETAY(RAD) MDDULO TWOP I 
THETAY(RAD) 

0 . 5A206A8E*0A A.5A20569 
0 . 15A2213E*02 2.855763A 
0 .563A020E 4 0A A. 2855981 



TDTAL CPU TIME = 8 MIN , 31.38 SEC 



TABLE 3.23 

LONG RANGE EIGENRAYS , PIECEWISE LINEAR SOLUTION 



RAY TRACING USING PIECEWISE LINEAR APPROXIMATION 
CASE ICHA 



NDATA = A 

G ( 1 ) = 0 . 16000E-0 1 l/SEC 

G ( 3 ) = 0 . 17000E-0 1 l/SEC 



NUMBER DF GRADIENTS = 3 

G ( 2) = -0.18889E-01 l/SEC 



ERAYS > T 




PRTALL 


= F 










Y0 = 100.0 M 




DEPTH = 


2000.0 M 


FREOC = 


250.0 HZ 






RNGSTP = 5.0 


M 


ANGSTP 


» 1.0 DEG 


1 HRZRNG 


» 50000.0 M 


YR = 250. 


0 M 


YERROR = 15.0 


M 






















EIGENRAYS 








RAY BETAO(DEG) 


DEPTH(M) 


TRVLT (SEC) 


BETAY(DEG) 


THETAT (RAD) 


THETAR(RAD) 


THETAY(RAD) 


MODULO TWDP 
















THETAY(RAD) 


1 A0. 000 


236.69 


52.A21823 


39.917 


0 .823AA0 1E*05 


0 . 3 36529 1E*05 


0 . A869109E*05 


2.689A315 


2 102.000 


261.00 


3A.AA0292 


102.535 


0.5A0?968E*05 


0 .5121057E*05 


0 . 2888 1 1 7E* 0A 


A. 1350528 


3 119.000 


238.33 


38.520812 


60.820 


0 . 6050835E*05 


0.A5790A0E*05 


0.1A71795E*05 


2.7282732 


A 122.000 


257. 7A 


39.73A666 


122.182 


0 . 62A 1507E*05 


0 . AA 3 9925E*05 


0. 1801581E*0S 


1.9203099 


TDTAL CPU TIME 


= A MIN 


, 27.61 SEC 
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RRY TRRCE-- PIECEWISE LINERR RPPROXIMRTION 
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a: « 
o 



( W)HlcDQ 

, IT 



Figure 3.15 Ray trace plot corresponding to Table 3.22 
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SSP 




Figure 3.16 



R ay trace plot 



corresponding to Table 



3.23 
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Table 3.23 shows simulation results using the piecewise 
linear method with a horizontal range (HRZRNG) of 50km. 
Parameters that affect execution time such as angle step size, 
depths, etc., match the simulation presented in Table 3.21 
using the ODE solver. Table 3.21 shows that the ODE solver 
used 24 minutes to complete a 200 meter range eigenvalue 
problem. In comparison. Table 3.23 shows that the piecewise 
linear method completed the 50km run in only 4 minutes. 

In summary, this chapter has shown that 

- the ODE solver performs accurate phase calculations 
with far fewer data pair samples when a smooth 
sound-speed profile is encountered; 

- increasing the range step size of the ODE solver to 
five meters reduces the CPU time cost without 
signficantly affecting the phase solution accuracies; 

- the piecewise linear method can run long-range 
simulations using relatively little CPU time; and 

- the piecewise linear method requires many sound-speed 
data pair samples for an accurate solution in a medium 
with a smooth, sinusoidal profile. 
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IV. CONCLUSIONS AND RECOMMENDATIONS 



The transfer function of the inhomogeneous ocean based on 
the WKB approximation requires solving a phase integral. Ray 
acoustics theory can provide solutions to this phase integral. 
The simulations performed show that the travel times 
calculated using the theories of ray acoustics can be used to 
solve the phase integral avoiding direct numerical 
integration. 

Two applications of ray acoustics produced computer 
simulation codes which 

- are capable of solving for the position, travel time 
and phase of a propagating ray, and 

- have very different advantages and costs. 

The first application was the piecewise linear approximation. 
Sound-speed versus depth data pairs sampled from the ocean 
medium were connected with constant gradient linear segments. 
Well-known, closed form equations form the mathematical model 
for sound propagation. The solutions are low cost (in terms 
of CPU time) , but many data samples are required for accurate 
phase solutions for arbitrary sound-speed profiles. 

In contrast, the Akima cubic spline/ODE solver method uses 
the medium samples to form a continuously variable sound-speed 
profile. Accurate phase calculations can be made with a 
minimum of sound-speed versus depth data pairs. The 
disadvantage in using the ODE solver is its exorbitant cost in 
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terms of CPU time. Simulations must be short range problems 
to limit computing costs. 

The ability to search for and identify eigenrays was 
developed for each of the two ray acoustics theory 
applications. Very distinctive characteristics were seen for 
phase calculations in terms of the number of data samples 
needed to assure accurate solutions. The identification of 
eigenrays is an easier task of position or depth computation. 
Here the distinctions tend to disappear with the two methods 
arriving at the same solutions for eigenray launch angles. 
The cost or CPU time required continues to strongly favor 
using the piecewise linear approach. 

These findings indicate that most propagation problems 
will require a piecewise linear approach for computational 
efficiency. Using the ODE solver would severely limit the 
range of computer simulations. 

Careful sampling of the ocean media is required to obtain 
accurate results from the piecewise linear method. Sufficient 
sound-speed versus depth data pairs must be used to accurately 
represent the sound-speed profile. 

This thesis developed the tools to quantify the strengths 
and weaknesses of two phase computation methods in a variety 
of media. With these findings in mind, future work 
recommendations are to 

- incorporate each phase computation technique as a 
module in the larger pulse propagation code; and 
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- run pulse propagation simulations to compare the 
received pulse shapes. 

If this comparison shows that the piecewise linear 
approximation produces a relatively undistorted received 
pulse, it is an efficient solution to the phase computation 
problem. 
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